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Abstract 



In this article we compare the mean-square stability properties of the 0-Maruyama and 
0-Milstein method that are used to solve stochastic differential equations. For the linear sta- 
' bility analysis, we propose an extension of the standard geometric Brownian motion as a test 

|^ , equation and consider a scalar linear test equation with several multiplicative noise terms. 

This test equation allows to begin investigating the influence of multi-dimensional noise on 
the stability behaviour of the methods while the analysis is still tractable. Our findings 
include: (i) the stability condition for the 0-Milstein method and thus, for some choices of 
6, the conditions on the step-size, are much more restrictive than those for the 0-Maruyama 
method; (ii) the precise stability region of the 6*-Milstein method explicitly depends on the 
, noise terms. Further, we investigate the effect of introducing partially implicitness in the 

' diffusion approximation terms of Milstein-type methods, thus obtaining the possibility to 

control the stability properties of these methods with a further method parameter a. Nu- 
merical examples illustrate the results and provide a comparison of the stability behaviour 
of the different methods. 
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^ '. 1 Introduction 

H | 

In recent years the area of numerical analysis of stochastic differential equations (SDEs) has 
expanded at a fast pace. This interest has been driven by different application areas, such as 
computational finance, neuroscience or electrical circuit engineering. A large part of research in 
stochastic numerics has been aimed towards the development and strong and weak convergence 
analysis of several classes of numerical methods. A further important issue for the investigation 
of numerical methods consists of examining methods for their ability to preserve qualitative 
features of the continuous system they are developed to approximate. A linear stability analysis 
is usually the first step of an analysis in this direction. For this the method of interest is applied 
to a scalar linear test equation and stability conditions on method parameters and step-size 
are derived and compared with the stability condition for the test equation. In deterministic 
numerical analysis the underlying idea for a linear stability analysis is based on the following line 
of reasoning: one linearises and centres a nonlinear ordinary differential equation x'(t) = f(t, x) 



* Support by the Leverhulme Trust through project "Stability issues for SDEs" (F/00 276/K) is gratefully 
acknowledged. 



1 



around an equilibrium, the resulting linear system x'(t) = Ax(t) (A the Jacobian of / evaluated 
at equilibrium) is then diagonalised and the system thus decoupled, justifying the use of the 
scalar test equation x'(t) = Xx(t), A G C, for the analysis. We refer to, for example, [3 Chapter 
IV. 2] for more detail on this procedure. 

In the stochastic case the same fundamental problem exists: we wish to preserve the qualitative 
behaviour of solutions of nonlinear stochastic differential equations following discretisation by a 
numerical method. Once again the starting point is a linear stability analysis. We can linearise 
and centre a nonlinear SDE around an equilibrium solution (see [10] for the corresponding the- 
ory), and this procedure yields a system of SDEs with an m-dimensional driving Wiener process 
of the form dX(t) = FX(t)dt + YT=i G r X(t) dW r (t), t > . Research on stability analysis for 
SDEs has focused on the scalar linear test equation dX(t) = XX(t)dt + fiiX(t)dWi(t), where 
its solution is called geometric Brownian motion. This corresponds to considering the linearised 
SDE system when it can be completely decoupled and with one driving Wiener process. Rel- 
evant references are given by, e.g., [3j [HJ [121 [20]. Some first explorations of a linear stability 
analysis for systems of SDEs have been made in [19[ 121] , and in [1] suitable linear test systems 
have been derived. 

The most common and well-known methods to treat SDEs numerically are the Euler-Maruyama 
approximation and the Milstein scheme developed in the last century. These methods and their 
drift-implicit counterparts have been studied for their stability behaviour, see for example [2"]lll[ 
[T2"| 120] . however always only for a one-dimensional noise term in the test equation. Our aim in 
this article is to investigate the influence of multi-dimensional noise on the mean-square stability 
behaviour of these numerical methods and to compare their stability behaviour. Assuming 
that the matrices F and G r can be diagonalised and the system above decoupled, we propose 
a scalar linear equation with an m-dimensional Wiener process, that is dX(t) = XX(t)dt + 
YlT=i^rX{t)dW r {t), as a suitable test equation. We consider the #-Maruyama and the 6- 
Milstein method and study the asymptotic mean-square stability properties of these methods. 
We find that the stability conditions for the 0-Milstein method are not only stronger than those 
for the #-Maruyama method, but they also become more restrictive when increasing the number 
of noise terms in the SDE. In particular we find it impossible to conclude A-stability for the 
^-Milstein method, in the sense that we can not define a value or bound Abound for the parameter 
9, such that the method applied to the test equation would be ^4-stable for all 6 > Abound for all 
m > 1 and all parameters \i r , r = 1, . . . , m. This indicates that analysing the stability behaviour 
of numerical methods using test equations with only a one-dimensional Wiener process may 
not provide sufficient insight into the properties of the methods for practical high-dimensional 
application problems. We then present a modification of the #-Milstein method, where we 
introduce partial implicitness involving a further parameter a in the hope to have a better 
control on the stability behaviour of the method. The stability analysis indicates that this is 
the case, although the stability condition still depends on the number of noise terms. Numerical 
experiments illustrate the influence that the choice of method and method parameters have for 
practical simulation tasks. 

In Section [2] we introduce the linear test equation, the 0-Maruyama method and the 0-Milstein 
scheme, we give a precise definition of asymptotic mean-square stability and state stability 
conditions for the test equation in the continuous case. Section [3] is devoted to the analysis of 
the stability properties of the methods in terms of the arising stochastic difference equations and 
we compare the stability regions of both types of methods in Section[H In Section[5]we introduce 
a new class of Milstein type methods and analyse their stability behaviour. We illustrate the 
theoretical findings by some numerical experiments in Section [6j 
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2 Preliminaries 



In this section we introduce the stochastic differential and difference equations, as well as the 
notions of stability, that we consider in this article. 

We will be concerned with asymptotic mean-square stability of the zero solution of a test equation 
with respect to perturbations in the initial data. As a test equation we employ the following 
scalar linear stochastic differential equation with multiplicative noise 

m 

dX(t) = XX(t)dt + f-i r X(t)dW r (t), t > t > 0, X(t ) = X , (1) 

r=l 

driven by an m-dimensional standard Wiener process W(t) = (W\(t), . . . , W m (t)) given on the 
probability space (f2, J-, P) with a filtration (J-t)t>t ■ By E we denote expectation with respect to 
P. We assume that the coefficients A and fi r , (r = 1, . . . , m) are complex- valued and without loss 
of generality suppose that the initial value Xo is non-random ( |16t Chapter 4]). The solution 
(X(t))t>t of (P) is a complex-valued stochastic process, where W is a vector-valued Wiener 
process in W 71 . Alternatively this complex- valued SDE can be written as 2-dimensional vector 
SDE in the real and imaginary parts. We denote the real and imaginary part of a complex 
number x by !EH(x) and 3(x), respectively, x denotes the complex conjugate and |x| stands for 
the absolute value of some x € C. 

For any non-zero initial value Equation ([1]) has a non-trivial path-wise unique strong solution 
which we denote by X(t;to, Xq) when we wish to emphasise its dependence on the initial data. 
For Xq = the equation obviously admits the zero solution, which is a steady state solution of 
the equation. 

Definition 2.1. (cf. ITty [76] / ) The zero solution of Equation ((TJ), X(t) = 0, is 

1. mean-square stable, if for each e > 0, there exists a 5 > such that 

E(\X(t;t ,X )\ 2 ) <e 

whenever t > to and \Xq\ < 5; 

2. asymptotically mean-square stable, if it is mean-square stable and if there exists a 5 > 
such that whenever \Xq\ < 6 

E(\X(t;t ,X )\ 2 ) ^0 for t ^ oo . 

The zero solution is called unstable if it is not stable in the mean-square sense. Definition 12.11 
is slightly more general than necessary in the present context, as for the simple linear equation 
given by ([1]) we can take S arbitrarily large and thus it does not play a significant role. 

Proposition 2.2. \1(A Ijgj/ The zero solution of Equation ([T]) is asymptotically mean-square 
stable if and only if 

1 m 

W + 2 J>r| 2<0 - ( 2 ) 

r=l 

We now discuss the stochastic difference equations that arise by applying the 0-Maruyama 
method and the 0-Milstein scheme to the scalar test equation ([T]). We consider numerical 
methods for computing approximations X/ « X(ti) of the solution of the test equation (JXJ) at 
discrete time points t% = ih (i = 0, 1, . . . ) with constant step-size h. 
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The #-Maruyama method applied to the test equation (UJ reads 

m 

X i+1 = X i + h(e\X i+1 + (l-6)\X i ) + VhY,VrX i Z r7 i, * = 0, 1, . . . , (3) 

r=l 

where we have replaced the Wiener increments j^'^ +h = W r (ti + h) — W r (ti) by the scaled 
random variables \fh£, r i. Here each {£ r ,i}ieN is ° ne of m independent sequences of mutually 
independent standard Gaussian random variables, i.e., each £ rj j is A/"(0, redistributed. 

The 0-Milstein method applied to the test equation (UJ is given by 

m 

x i+l = x i + h(exx i+1 + (i-e)xXi) + VhJ2fjt r x i ^ i 

r=l 

^ m ^ m 

+ 2 h ^2^r X i (£r,i ~ !) + 2 k ^2 ^n^Xi^^r^i, % = 0, 1, . . . . (4) 

r=l n> r 2 =1 

Here we have additionally replaced the multiple Wiener integrals 

I^+ h := / 1 / dW r2 (u)dW ri (s) 
Ju Ju 

as follows: (i) If r\ = T2 the multiple Wiener integral Ir l f l+h can be replaced by ^(£ri ~~ -0 f° r 
r = 1, . . . , m, and (ii) if r\ ^ T2 we use the identity I ri ,r 2 + Ir 2 ,r 1 = Ir 1 Ir 2 to obtain 

fl ri /J lr2 I rir2 + fl r2 fl ri I r2 ri ll ri [l T2 (l ri T2 -\- Ir2,ri) l^n t^ra^Ti Ira /^riMr2^£r*i,i£r2,i ■ 

For = the method (|3j) reduces to the (forward) Euler-Maruyama scheme, which is explicit 
in the drift as well as in the diffusion part. For 9 > the methods © and ([U are drift-implicit. 
For 9 = 1/2 the scheme (|3|) is known as stochastic trapezoidal rule and for 9 = 1 we obtain the 
backward Euler-Maruyama method. 

The two last terms in the #-Milstein method represent a higher order approximation of 
the diffusion part in Equation (pQ) in the sense of mean-square convergence. We call a method 
mean-square convergent with order 7 (7 > 0) if the global error, X(ti) — X$, satisfies 

max \\X(ti) - Xi\\ L2 := max (E\X(U) - Xi\ 2 ) 1/2 < CK 1 as h->0, 
i=0,l,... i=0,l,... 

with a positive error constant C, which is independent of the step-size h. It is well-known that 
in the case of multiplicative noise the 0-Milstein method @ is mean-square convergent of order 
7 = 1, whereas the 0-Maruyama method ([3]) is mean-square convergent of order 7 = 1/2. 

Obviously the stochastic difference equations ([3]) and (|4|) admit the zero solution for the initial 
value Xq = 0, which is a steady state solution as well. For any non-zero initial value Xq, the 
equations have a unique solution provided (1 — hX9) 7^ 0. We write Xi(to,Xo) again when we 
want to emphasise that the solution of the difference equations depends on the initial data. 

Definition 2.3. The zero solution of the difference equations ([3]) and ([4]) is 

1. mean-square stable, if for each e > 0, there exists a 5 > such that 

E(|X i (t ,X )| 2 )<e 

whenever i > and \Xq\ < 5; 
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2. asymptotically mean-square stable, if it is mean-square stable and if there exists a 5 > 
such that whenever \Xq\ < 5 

E(\Xi(t ,X )\ 2 ) -+0 for i^oo. (5) 

Stability conditions will now involve the coefficients A, /j, r (r = 1, . . . , m) of the test equation 
dH), as well as the method parameter 9 and the step-size h. 

It will be useful to describe the stability region of a stochastic differential or difference equation. 
We follow the presentation in [12] and consider sets of parameters S'sde, Se-Marf^ h), Sg_Mii(9, h) 
for which the zero solutions of the continuous and the discrete equations are asymptotically 
stable, that is 

-S'sde := {A,/ii, . . . , y, m € C : ©holds}, 
Se-M&r(8, h) := {A, fii, . . . , fi m £ C : holds for solutions of p}}, 
S'e-Mii(^) ^) := {A, /tii, ... , /i m € C : (jSJ) holds for solutions of (jl}}. 

Further, we consider the extension of the deterministic notion of A-stability [£l [15] to the mean- 
square analysis setting and say that a numerical method is A- stable in mean-square |11[ 112]. if 
whenever the zero solution of ([I]) is asymptotically mean-square stable, then the same is true 
for the zero solution of the method for any step-size h > 0. Using the above definitions of the 
stability regions, we call the #-Maruyama or #-Milstein method A-stable in mean-square if for 
all h > we have S'sde C S d _ Ma ,r(0, h) or S'sde Q S _ M ii(#, h). 



3 Stability analysis for the numerical methods 

For the stochastic difference equations and (j3J) we now derive stability conditions in depen- 
dence on the method parameter 9 and the applied step-size h, and compare these conditions 
with those for the continuous problem given in Lemma [ 



3.1 The 6>-Maruyama method 

We start by rearranging the stochastic difference equation ([3j) into the following one-step recur- 
rence equation. Then by squaring and taking expectations we obtain a recurrence equation for 
the second moments E|Aj| 2 . As mentioned before, we have to assume that (1 — QhX) ^ to 
guarantee the existence of a unique solution of the recurrence equations. 

Rearranging Equation © yields the recurrence 



X i+l = I a + ^2 b r £r,i J X i 
\ r=1 / 



0,1,..., (6) 



where 

h\ Vhfi r 

a:=1 + T^9Tx and br:= T^eT\- (7) 

Then the recurrence equation for the second moment E|Aj| 2 , using E(£ rj j) = 0, E(£ 2 J = 1, and 
the complex conjugate equation of ([6]), reads 

(m m \ 

(a + ^ b r £ r ,i) (S + l r ) E l^| 2 
r=l r=l ' 

m \ 

a| 2 + £|6 r | 2 )E|^| 2 . (8) 

r=\ ' 
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From ([5]) we can immediately read off necessary and sufficient stability conditions in terms of 
the above parameters. 



Lemma 3.1. The zero solution of the recurrence equation ©, representing a one-step Maruyama- 
type method applied to the test equation (pQ), is asymptotically mean-square stable if and only 
if 



|a| 2 + ^|6 r | 2 < 1. (9) 



r=l 

Now rewriting the stability conditions Q in terms of the parameters A and fj, r (r = 1, . . . , m) 
in the test equation ([1]), the method parameter 9 and the step-size h, using (J7J), we obtain the 
following result. 

Corollary 3.2. The zero solution of the stochastic difference equation given by the 6-Maruyama 
method (|3j) applied to the scalar linear test equation §1$ is asymptotically mean-square stable if 
and only if 



1 m 1 

^(A)+2El^! 2 + 2 /l(1 " 20)|A|2 < °- (10) 

r=l 

We note that the first two terms in the left-hand sides of ()10p are equal to the left-hand side 
of ([2]), that is they correspond to the stability condition for the continuous problem. 

Now comparing the stability condition for the continuous problem to that of the discrete problem, 
we immediately obtain from (jlOp an extension of the result [121 Thm. 4.1] to the case of ([I]) 
driven by a multi-dimensional Wiener process. 

Corollary 3.3. For all h > it holds that 

So-mrfrh) C Ssde for 0<9<l/2, 
Se-M&r{9,h) = S SDE for 6 = 1/2, 
Se-M,r(0,h) d Ssde for 0>l/2. 

In particular, for 9 > 1/2 the 9-Maruyama method is A-stable in mean-square. For < 9 < 1/2 
and (A, m, . . . ,/i m ) S Ssde? the stability condition ([TO]) for the 9-Maruyama method is satisfied 
if and only if 

-2ffl(A) + iEr=iKI 2 ) m) 

< (1 - 2(9)|A| 2 ' [ ' 

3.2 The 6>-Milstein method 

We now turn to the #-Milstein method and first follow the same steps as in the previous section 
to deduce a recurrence equation for the second moments of the solution of (j4]) and read off the 
corresponding stability conditions. Rearranging the difference equation (j3|), we obtain 



(mm \ 
Oj ~\~ ^ ^ b r + ^ ^ Cj- i r2 £ri,i£r2,J ) 
r=l n,r - 2=l 



(12) 



in 



\fh[i r 2 h(j, ri fj, r2 
where a := a - 2_^ c r,r, K := _ c Tl . X2 = - _ — , (13) 



r=l 



G 



and o is given in ([7]). Note that rewriting the parameter o in terms of a provides a convenient 
way to compare the stability conditions for the #-Maruyama and #-Milstein method. 

With analogous calculations as in the previous section and by additionally using E(£^) = and 
E(v i) = 3, we find the recurrence for the second moments of the 0-Milstein method as 

E|^i+i| 2 



( m in m in \ 

( Q ~l~ ^ 1 br £r,i ~l~ ^ ~] Cri,r2 £ri,i£r2,i\ ( a ^""^ Cr,t ~l~ ^ ^ c ri,r2 £ri,j£r2,i^ j E|-^Q| 
r=l ri,r2=l r=l ri,r'2=l 



ri,r2=l r=l n.j»"2 = 

m mm mm 

oa + a ^ c r>r + a ^ cv.r + ^ f>A + ( ^ Cr ' r ) ( 5^ 5r > r 

r=l r=l r=l r=l r=l 

m mm 



in in in \ 

^ ^ Cr,r-Cr,r ~\~ ^ ^ Cri,r'2 ■ ^ ^ *-Vi,r 2 j E|-Xj| 

m m m m mm 

+ ^ ^ c rr J ( a + ^ ^ Cy r J ~\~ ^ ^ -\- 1 ^ ^ c r)r c rjr - + ^ ^ c ri i? - 2 • ^ ^ c ri r2 ] E | -X"; 



r 1 ^r 2 r l^ r 2 

m mm 

a 

r=l r=l r=l r=l r 1 ,r 2 =l r 1 .r 2 = i 

r 1 ^r 2 r li ZT 2 



2 

% 



(m m m m \ 

| CL | -|- ^ ^ 1 6 r | -|- 2 ^ ^ | C r)J - 1 + ^ ^ Cr\ ,r2 ' ^ ^ ^Vi ,r2 ) 
r=l r=l r 1 ,r 2 =l r 1 ,r 2 = l / 
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(14) 



Again we can read off from (|14[) necessary and sufficient stability conditions in terms of the 
parameters a, b r and c Tl)T2 . 

Lemma 3.4. The zero solution of the recurrence equation (|12p given by the 6-Milstein method 
applied to the test equation ([I]), is asymptotically mean-square stable if and only if 

m m 

|a| 2 + X^| 2 + 2 XX^ 2 + l c ^ m l 2 < *> ( 15 ) 

r=l r=l 

where c sum = ^2 rijr2 =i jri jL r2 c rijr2 . 

Remark 3.5. When applied to the linear scalar stochastic differential equation ([I]), other vari- 
ants of one-step Maruyama-type methods and one-step Milstein-type methods can quite often be 
rearranged into the recurrence equations §6§ and (|12|) . respectively, with appropriate definitions 
of the parameters a, b r and c ri)T2 . Then Lemmata \3.1\ and \3.4\ can be applied to obtain stability 
conditions interpreted with these definitions of a, b r and c rijJ . 2 . 

The next corollary follows by rewriting the stability condition f)15j) in terms of the original 
parameters using (fl~3j) . 



Corollary 3.6. The zero solution of the stochastic difference equation given by the 6-Milstein 
method (jlj) applied to the scalar linear test equation (TjQ) is asymptotically mean-square stable if 
and only if 

^(A) + -^|^| 2 + ^Ml-2^)|A| 2 + -/ i ^|/i r | 4 + -/ i I /Wr 2 < 0. (16) 



2^' r " 2 v " 1 4 

r=l r=l 
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Again the first two terms in the left-hand side of (|16|) are equal to the left-hand side of ^) . How- 
ever, when comparing the above condition with the stability condition (|10p for the #-Maruyama 
method, we see that the left hand side of (|16|) contains two additional terms, which are always 
non-negative and depend on the noise intensities fj, r (r = 1, . . . ,m) but are independent of the 
parameter 9. Thus the precise stability region Sq_mu(0, h) of the 0-Milstein method depends on 
the noise intensities \i r (r = 1, . . . , m). Suppose our aim is to use the method parameter 9 to 
determine the optimal stability region of the #-Milstein method for a given set of parameters 
in ([1]), that is we aim to find a value of 9 such that the sum of the third to the fifth term in (|16p 
vanish for any step-size h. Then, assuming that (j2j) holds, we obtain for every set of parameters 
Hi, ... , Hm a different value for that optimal # Mlj ... iAim . We define 



'/«.,. ..,/Xti 



1 4 z 

- + V lM + l/Xsum| (17) 

2 ^4|A| 2 ^ 8|A| 2 ' 1 ; 



where 



(18) 

For this optimal On^ ^ we immediately obtain from the stability condition (|16|) for the 0- 
Milstein method the following corollary. 

Corollary 3.7. For all h > it holds that 



S e .Mn(0,h) C Ssde for 0<9< 
Se-mi(Q,h) = S S de for 6 
Se-Mii(0,h) D Ssx>£ /or > 6^ u 



'/ii,...,ti„ 



Assuming that (A,/ii, . . . , € Ssde? /or < 9 < 0u\ ...,n m i the stability condition f)16|) /or i/ie 
9-Milstein method is satisfied if and only if 



-2(*(a)+±e; 



;i-2^ia|2 + ie"LiKi 4 + ii// 



simi 



In particular, condition (JT9J) implies that for 9 > \ and (1 — 2#)|A| 2 + ^ E^=i l/ i r| 4 + \\^sum\ 2 < 0, 
t/ie stability condition (|16|) /or t/ie 9-Milstein method is satisfied for any step-size h > . 
In i/ie case £/ia£ (A, H\, . . . ,pL m ) Ssde, and 9 > \ and (1-26>)|A| 2 + \ Y^=\ l/^l 4 : + \\v-sum\ 2 > 0, 
then the zero solution of the 9-Milstein scheme is also unstable only if condition (I19p is imposed 
on the step-size h. 



Intrinsic in the concept of ^-stability is the idea that the property of A-stability of a method 
holds for the whole class of differential equations considered as test equations. In the setting of 
this article this implies that we would need to find a 0bound> a value of or a bound on 0^,...,^, 
which is independent of A, fj,±, . . . , fi m , such that Sq_mu{9, h) 2 Ssde for 9 > 6*b OU nd- In the case 
of the #-Maruyama method the corresponding value of Abound is \, see Corollary 13.31 
However, considering the simple case of multi-dimensional noise terms having equal noise in- 
tensities, we find that, at best, we can find an upper bound Abound independent of the given 
parameters A, . . . ,fJ, m , but depending on the number of noise sources. To see this, let fii = 
fi2 = • • • = fJ-m- Then, using the squared stability inequality (J2J) in the form ||a*i| 4 /|A| 2 < 
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m 

and Yl Mi/^i = (m 2 — m) in (fT7|) yields 

n.,r2=l, r\^r% 



< ^ + - + ^(m-l) 2 , (20) 
2 m 2 

Thus for example, for one, two and three noise sources (m = 1, 2, 3), we can find a Abound a s 
| , | and -jr , respectively. In the one-dimensional noise case this corresponds to the result in jilt 
Cor. 2.2]. But in general there exists no upper bound Abound for 0^ jfim such that the 0-Milstein 
method for any 9 > Abound is ^4-stable for the whole class of equations ([T|) with arbitrary many 
noise terms. In the general case of the SDE ([T]) it is possible, essentially using the Cauchy- 
Schwarz inequality several times, to derive an upper bound Abound for ^i,...,^ m which also only 
depends on m, but as it involves m 3 it becomes pointless for large m. Thus, it appears that 
the stability condition (|16p for the 0-Milstein method becomes very restrictive for an increasing 
number of noise terms. 



4 Comparisons of the stability regions 

In this section we aim to illustrate the results of the previous sections by visually comparing the 
stability regions Se-Mar{0-,h) and S$.Mii{0,h). As there are too many parameters in the system 
to do so in a two-dimensional plot, we only consider the case of the test SDE (pQ) with multi- 
dimensional noise where all the terms have the same noise intensity, i.e., /Ui = [ii = ■ ■ ■ = fi m , 
and real- valued coefficients, i.e., we have A, fj,\ € R. We essentially follow [12] regarding the 
scaling of the parameters in the plots. Thus we set 

x := hX and y := hm/if . 

The stability conditions fl2J), (fTUj) and (fT6|) become 

x+\y < 0, (21) 
x + l y+ ] -(l-26)x 2 < 0, (22) 
x + ^y + ^(l-29)x 2 + ~y 2 + ^y 2 (m-l) 2 ^ < 0. (23) 

Note that for m = 1 and m = 2 condition ()23p yields the same stability region. 
To interpret the figures below, observe that given a test equation with parameter values A and fi\ 
and m = 1, the point (x, y) = (A, n\) corresponds to the choice of step-size h = 1. Then varying 
the step-size h corresponds to moving along the ray that connects (A, fj, 2 ) with the origin, where 
going on this ray in the direction of the origin corresponds to decreasing the step-size h. For 
m > 1 the scaling is appropriately adapted. 

Figure [T] shows the mean-square stability regions of the zero solutions of the SDE (white area 
with a dashed border), the 0-Maruyama approximation (light-grey area) and the #-Milstein 
approximation (light dark area to dark-grey area) for different values of the method parameter 
9 and illustrates how the stability region of the #-Milstein method decreases with m and compares 
to that of the #-Maruyama method. In particular, one can observe that the mean-square stability 
regions for the #-Milstein method are always smaller than the stability region of the 0-Maruyama 
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method. For the 0-Maruyama method the figure illustrates the property of ^4-stability of the 
method for 9 > ^ for all m > 1, where as for the 0-Milstein method we can not conclude 
^4-stability for a particular 9 and all m > 1. 



U = 6 = 0.25 e = o.5 




-10 -5 5 10 -10 -5 5 10 -10 -5 5 10 



Figure 1: SDE |JT]) with n\ = fa = • • • = Mm and x = hX and y = hmfj,f: Mean-square stability regions 
for the zero solutions of the SDE (white area with dashed border), the 0-Maruyama method (light-grey 
area) and the Milstein method for m = 2, 3, 4, 5 (light dark-grey area to dark dark-grey areas). 

In [12^ Section 4] and Section 3] the stability regions of the #-Maruyama method and the 
0-Milstein method have been plotted separately for the scalar SDE (prj and real- valued X,fJ,\. 
To emphasise the different stability properties of both methods even in the case of m = 1, 
we provide in Figure [2] a plot of the stability regions of both methods together for the same 
setting as in [1 2 1. [TT] and several values of the method parameter 9. In these cases the stability 
regions of the 0-Maruyama method and the #-Milstein method coincide with that of the SDE if 
9 > 1/2 and 9 > 3/2, respectively. When using the Euler-Maruyama method and the standard 
Milstein method, that is taking 9 = 0, the methods that are most often used, then it appears 
from Figure [2] that both stability regions are quite small but of similar size. This is already 
mentioned in [21 p. 114]. However, when taking mean-square accuracy into account and aiming 
to reduce numerical costs by using the Milstein method with a larger step-size, one might easily 
be prevented from doing so by the more restrictive stability condition of the Milstein method. 

5 The effect of introducing partial implicitness in the discreti- 
sation of the diffusion term 

A brief look at the deterministic case, that is \i T = (r = l,...,m) in ([T|), (jiop . etc., 
reminds us that it is by making the approximation implicit and the introduction of the method 
parameter 9 that one can control the stability properties of the method by suitably choosing 
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-10 -5 5 10 -10 -5 5 10 -10 -5 5 10 



Figure 2: SDE (p} with m = 1 and x — hX and y = hfij: Mean-square stability regions for the zero 
solutions of the SDE (white area with dashed border) , the 6*-Maruyama method (light-grey area) and the 
0-Milstein method (dark- grey area). 

9. In this sense it would be useful to develop Milstein-type methods that incorporate implicit 
higher order approximations of the diffusion term to have a method parameter for this purpose. 

In general, a straightforward introduction of implicitness into approximations of the diffusion 
term results in the numerical solution becoming unbounded with positive probability, see |18| 
Chap 1.3.4] and [8] for methods avoiding this problem. In this article our aim is to highlight the 
effect that an implicit discretisation of the diffusion term can have on the stability properties of 
the 0-Milstein method and thus we take advantage of the simple structure of the test equation ([I]). 
The latter results in the fact that in each term in the second last sum in the 0-Milstein method ([4]) 
the double Wiener integral J f s dW r (u) dW r (s) can be replaced by 5/1(^,1 ~~ 1)> for r = 

1, . . . , m. Then the second last term in the method ([!]) can be written as \ hp^Xi £^ — | h/J^Xi 
for each r and we introduce an implicit approximation with an additional positive method 
parameter a only in the latter term, which does not contain the random variable £ r j. Thus we 
propose, for each r, to use 

- hf? r (Xi ££i - (<rX i+ i + (1 - a)Xi)^ instead of - hfij (Xi - X^j. 

We emphasise that this represents a (partial) implicit approximation of the diffusion term in 
contrast to well-known approaches to implicit approximations of the drift term. 
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Applied to the scalar linear test-equation (pTJ) the #-<r-Milstein method then takes the form 

m 

X l+1 = X i + h(0\X i+1 + (l-e)\X i )+VhJ2VrXitr,i 

r+l 

rn m 
+ 2 h J2Vr( X i£,i ~ i aX i+l + i 1 ~ a ) X i)) + 2 k Yl ^rx^r 2 X^ ri 4 r2>i , (24) 

r = l r 1 ,r 2 = l 

i = 0, 1, . . . , with the random variables £ rj j defined as in Section [2j 

Remark 5.1. The choice of a = 9 yields the 8-Milstein approximation of the Stratonovich-SDE 
dX(t) = (XX (t) - \ £^=1 tfX(t))dt + £™i ^X{t) o dW r (t). 

For the stability analysis we follow the same procedure as in Section and first rewrite (|24|) as 
a one-step recurrence 



X i+1 = — [ a, + ^ ^ b r ^ rj i + ^ ^ C rijr - 2 £ri,i£r2,i ) -^i: 
\ r=l ri,r2=l 



(25) 



ri,r2= 

where the parameters a, b r ,c ri)T2 and (i are given by 



^ m 

a := 1 + (1 — #)/iA — 2^(1 ~~ CT ) ^ ' ^ r := ^^Mr > 

r=l 

1 1 m 

-/i/i ri/ u r2 , d := I — 9h\ + -ah'^2 /j% . 



c ri,r2 •" 2' 



r=l 



We define a := a + X^^Li c r,r and assume that (1 — 6h\ + ^cr/i^^Li A* 2 ) / to guarantee the 
existence of a solution to Equation (125p . Now squaring and taking the expectation yields 

E|X m | 2 = p|2^|a| 2 + ^|6 r .| 2 + 2^|c r , J .| 2 + |c^ m | 2 J E|X;| 2 , (26) 
where c sum = ^2,r\.T 2 =\,T\^r 2 c n,r 2 - 

Hence the zero-solution of the stochastic difference equation (|25p is asymptotically mean-square 
stable if and only if the factor on the right hand side of (|26p is less than 1. Rewriting this 
condition in terms of A, /i r , h, 9, a and rearranging, we obtain the following result: 

Lemma 5.2. The zero solution of the stochastic difference equation given by the 9-o-Milstein 
method (|24p applied to the scalar linear test equation ([1]) is asymptotically mean-square stable if 
and only if 

K ( A ) + 2 H 2 + 2 k[l ~ 20)|A|2 + 4 h Y^ + 8 h ^™m\ 2 + -ahJ2^l) < , (27) 

r=l r=l r=l 

w/iere a#am fi sum : = En,r 2 =i,r 1 ^r 2 AVi/V 2 - 

We note that the first terms in the left-hand side of (j27|) are equal to the left-hand side of (|16p . 
The additional term 2 a ^^2r=i 5?(A/x 2 ) is negative for 3?(A) < 0, i.e., when the stability condi- 
tion (|2|) for the test equation ([1]) is satisfied. Thus, the stability condition (|27p in Lemma 15.21 is 
less restrictive than the condition (1161) for the #-Milstein method. 



Denoting by Se_ a _Mii(9, cr, h) the stability region of the #-<7-Milstein method, analogous consid- 
erations as for Corollary 13.71 yield the next result. 
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Corollary 5.3. Define 



1 | Mr I , l/Wt|_ , ^ (^Mr 

" 2 ^4|A| 2 8|A| 2 2 |A 



For all h > it ZioZds £/ja£ 



[ll,...,fj,m ) 



5e. (T -Mu(f,o' ) ^) C SsDE /or 0<6> + <7< 
Se-<r-mi(0,<7,h) = S S de for 9 + a = 6» w ,... !(Um , 
So- a -Ma{Q,0,h) D S S de for 9 + a > 9 ' ill ,..., IMn . 

Then, assuming that (A, Mi, . . . , fi m ) G .Ssde, /or < + cr < f^,...,^ i/ie zero solution of the 
9-a-Milstein method is asymptotically mean-square stable if and only if 

. < -2(W(A) + | E"=i K| 2 + |o E^i K(A^)) f . 

(l-20)|A| 2 + ±£™ J/^ + ll/wl 2 ' 

In particular, condition (128p i/ien implies that for 6 > ^ and (1 — 2#)|A| 2 + | Sr^i lMr| 4 + 
jlMsiiml 2 < 0, i/te stability condition ([27]) /or i/ie 6-Milstein method is satisfied for any step-size 
h>0. 

As in Section U] we illustrate the stability regions of the #-o~-Milstein method applied to the test 
SDE ([!]) with a single noise (m = 1) in terms of real-valued coefficients A, /ii G K- In this case 
we can find a bound on 9u 1 ,...,u m as 3/2. Again, we set x = hX and y = /i//?; an d the stability 
inequality (f27l) in terms of x and y reads 

x + + 1(1 - 20)x 2 + ly 2 + ^axy < . 

Figures [3] - compare the stability regions for the SDE ([T]) with m = 1, the #-Maruyama 
method and the #-o~-Milstein method. In each Figure we have fixed the parameter 9 and only 
the parameter a varies. For a = the plots would correspond to those in Figure [2J We can see 
that the stability region of the #-o~-Milstein method (dark-grey area) is becoming larger when 
increasing the parameter a. Moreover, the stability region coincides with the stability region of 
the test equation ([TJ (white area) if# + o~>3/2. 




-10 -5 5 10 -10 -5 5 10 -10 -5 5 10 



Figure 3: Mean-square stability regions for the zero solutions of the SDE (white area), the 0-Maruyama 
method (light-grey area) and the 6>-er-Milstcin method (dark-grey area). Here 9 = and a = 0.5, 1, 1.5. 
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o = a = 0.5 a = 1 




-10 -5 5 10 -10 -5 5 10 -10 -5 5 10 



Figure 4: Mean-square stability regions for the zero solutions of the SDE (white area and dashed line) , 
the 0-Maruyama method (light-grey area) and the 9-a- Milstein method (dark-grey area). Here 9 = 0.5 
and a = 0,0.5, 1. 



o = o = 0.5 o = 1 




-10 -5 5 10 -10 -5 5 10 -10 -5 5 10 



Figure 5: Mean-square stability regions for the zero solutions of the SDE (dashed line), the (9-Maruyama 
method (light-grey area) and the ^-c-Milstein method (dark-grey area). Here 9=1 and a = 0, 0.5, 1. 



e = o e = 0.25 e = 0.5 




-10 -5 5 10 -10 -5 5 10 -10 -5 5 10 



B=1 8 = 1.5 6 = 2 




Figure 6: SDE ((T|) with /zi = /12 = • • • = [i m and x = hX and y = hm^i\: Mean-square stability regions 
for the zero solutions of the SDE (white area with dashed border), the 0-Maruyama method (light-grey 
area) and the 9 — cr-Milstein method for a = 1 and m = 2, 3, 4, 5 (light dark-grey area to dark dark-grey 
areas). 
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Further, we consider again the SDE ([T]) with identical noise intensities fix = [i<i = . . . = /i m and 
plot in Figure [6] the stability regions for m = 2, . . . , 5 (again the condition below is the same 
for m = 1 and m = 2) and a = 1 corresponding to the scaling x = hX and y = hmfif, where 
condition (f2"T|) now reads 

x + rl/ + r(l - 20 x 2 + 7 -y 2 + ^ m - 1 + -axy < . 
2 2 4 m 8 2 

Again the stability region of the #-<j-Milstein method decreases with growing m, as in the case 
of the 0-Milstein method, which prevents to conclude ^4-stability for a particular value of 9 + a 
and all m > 1. However, the introduction of the (partial) implicitness in the diffusion term 
clearly provides an improvement of the stability behaviour for the #-<7-Milstein method. 

6 Numerical experiments 

In this section we aim to illustrate the effect that the choice of method and the choice of the 
method parameters 9 and a have on practical simulation runs. We consider the test equation 
([1]) with the following parameters: m = 3, A = —2, fj,i = 1,^2 = — 1,^3 = 1, = 0.1, and 
integrate for t £ [0, 10]. The equation has the explicit solution X(t) = 0.1 • exp((— 2 — |)t + 
W 1 (t)-W 2 (t) + W 3 (t)). 

The following four figures show numerical simulation studies performed with the #-Maruyama 
method (Figure [7]), the #-Milstein method (Figure [8]), the #-<7-Milstein method with a = 1 
(Figure [9]) and the #-cr-Milstein method with a = 1.5 (Figure [10|) . all simulations done with 
a fixed step-size h = 1 over the interval [0, 10] with gridpoints U = ih, i = 0, . . . , 10. The 
parameter 6 varies as 9 = 0, 0.5, 1, 1.5. 

Remark 6.1. The analytical results in the previous sections also suggest choices of step-sizes, 
such that the zero solution of any one of the numerical methods for a fixed set of parameters 
is asymptotically mean-square stable and obviously it is possible to illustrate this by performing 
numerical experiments with fixed parameter sets and varying the step-size h. However, the focus 
of this article is rather on comparing qualitative properties of methods than on how individual 
methods behave for various step-sizes. 

Each of the four figures below consists of three pictures. The underlying idea of a stability 
analysis of numerical methods is to provide guidance for choosing method parameters and a 
step-size such that the numerical solution represents a good approximation of the true solution. 
(Note that convergence of a method only guarantees this in the limit for h — > 0.) Thus, we plot 
in the left picture of each figure the mean-square error e between the exact solution X(ti) and 
numerical solution X{ for the corresponding method and choice of method parameters. This 
allows to compare the impact of the choice of method on actually computing solutions of SDEs. 
However, the quantity of interest in the stability definitions 12.11 and 12.31 are the second moments 
of the analytical and numerical solutions. Therefore the other two pictures present plots of the 
second moments of the solutions of ([T]) , ([3]) , @ and ()24|) , in the middle picture estimated from 
the numerical simulations of ([3]), (j4]) and (|24j) . in the right picture computed form the analytical 
solutions. This allows to compare the effect of the choice of method and parameters on the 
behaviour of the second moments analytically and also, how well this is approximated by the 
numerical methods. 

In detail, the figures present the following quantities, where (X(ti,uij))i=o t ,,, t N and {Xi{ujj))i = \ : ,,, t N 
denote the values on grid points of a trajectory of the explicit solution of (pQ) with the above 
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parameters and of the numerical trajectory produced by one of the methods, respectively: 



(M \ 1 / z 

— ^(Xfauj) - Xiiuoj)) 2 « (E\X(ti) - Xtf) 1 ' 2 , 

1 - 

E ( X ") = mE^K) foralH, (29) 
i=i 

E(X 2 (ti)) = X 2 -exp j(2A + ^ / u 2 )-i i l foralH! (3 ) 



r=l 



E(X/) = s l ■ X^ for all i . (31) 

The expression (|29p for E(AC 2 ) represents an estimator for the second moment of the numeri- 
cal approximation, the expression flU for E(X 2 (^)) is the solution of the deterministic ODE 
E(X 2 (t))' = (2A + YT=i f4) E ( x2 ( t )) ( see P roof of Lemma EL_1 at discrete time-points and the 
last expression for E(Xf) is the result of applying the recurrences ([8]), (|14p and (|26p where s 
denotes the factor in front of E(Xf) in the right-hand side of the corresponding recurrence equa- 
tion. The number of trajectories computed for the above quantities is M = 100000 for each 
simulation. 

Figures [7] to [10] illustrate the behaviour of the methods for the above set of parameters, in 
particular the same fixed step-size h = 1 and different choices of 6 and a. The #-Maruyama 
method provides reliable approximations for 9 > 0.5, but not for 6 = 0, the #-Milstein method 
in same setting does not provide reliable approximations for 9 = and 6 = 0.5. In fact, the 
numerical solutions for 9 = 0.5 in this setting diverge. We can observe the improvement in the 
stability behaviour when using the #-cr-Mhstein method introduced in the previous section. In 
contrast to the #-Milstein method this method produces reliable approximations for the choice 
of 9 = 0.5 when also setting a = 1. Further, for the choice of a = 1.5 the #-cr-Milstein 
approximation behaves satisfactorily for all values of 9, thus the implicit term in the diffusion 
approximation is effectively stabilising the method. 



0.01 
0.009 
0.008 
0.007 
0.O06 
0.O05 
0.004 
0.003 
0.O02 
0.001 



MS-error over 1 00000 sample paths 



8=0 
8=0.5 




estimated second moment over 100000 sample paths 



8=0 
8=0.5 



analytic second moment 



2 4 6 




Figure 7: MS-error, estimated and analytic second moments using the fl-Maruyama method with step- 
size h = 1. 
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Figure 8: MS-error, estimated and analytic second moments using the 0-Milstein method with step-size 
h = 1. 
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Figure 9: MS-error, estimated and analytic second moments using the 0-tr-Milstein method with step-size 
h = 1 and a = 1. 
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Figure 10: MS-error, estimated and analytic second moments using the f/'-er-Milstem method with step- 
size h — 1 and a = 1.5. 
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7 Conclusions 



A linear stability analysis has been performed for the #-Maruyama method and the #-Milstein 
method, using a linear test equation with several multiplicative noise terms. We have obtained 
stability conditions guaranteeing asymptotic mean-square stability of the zero solution of the 
stochastic difference equations resulting from both types of methods applied to the test SDE. 
Comparing the stability conditions (jlOp and (|16p for the #-Maruyama and the #-Milstein method, 
respectively, it is quite obvious that the latter is more restrictive than the former, due to the 
fourth and fifth term in (|16p . i.e. \hY^ = \ |/v| 4 + \h\fJL sum \ 2 , which is always non-negative. In 
particular, this term is the more restrictive, the larger the parameters in the diffusion term in (TjQ) 
are. Now, in the case that these parameters are small, that is the so-called small noise case, it 
is well known from the corresponding mean-square convergence analysis that Maruyama-type 
methods provide sufficiently accurate methods for practicable choices of step-sizes, see [6l fTT| [22] . 
Milstein-type methods include higher order approximations of the diffusion term with the aim 
that they are accurate methods for more practicable choices of step-sizes just in the case that 
the diffusion term is large. In other words, when dealing with SDEs with larger diffusion 
terms, numerical efficiency considerations would suggest using a Milstein-type method with a 
larger step-size rather than a Maruyama-type method with a small step-size to obtain numerical 
approximations of the solution of an SDE with a similar accuracy. However, Corollary 13.61 
indicates that in this case there may be a trade-off and one is faced with restrictions on the 
step-size for the Milstein-type method due to stability reasons. Further, we have shown that 
the precise stability region of the #-Milstein method depends on the number and magnitude 
of the noise terms, whereas the stability region of the 0-Maruyama method is independent of 
them. In particular, it is not possible to define a value Abound such that the 0-Milstein method 
can be called j4-stable for all 9 > Abound for the class of SDEs (HJ for an arbitrary number of 
driving Wiener processes. For the #-Maruyama method Abound = \, as in the deterministic case 
and the case of the test equation ([1]) with m = 1. We provide a modified Milstein-type method 
with a partially implicit diffusion approximation and demonstrate that the resulting stability 
behaviour can be controlled more favourably. The results highlight that it is necessary to include 
multi-dimensional noise into test equations and to study their effects on the practical behaviour 
of the methods. 
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